module glcBehaviorMod

  !-----------------------------------------------------------------------
  ! !DESCRIPTION:
  ! Determines a number of aspects of the behavior of glacier_mec classes in each grid
  ! cell.
  !
  ! !USES:
#include "shr_assert.h"
  use shr_kind_mod   , only : r8 => shr_kind_r8
  use shr_log_mod    , only : errMsg => shr_log_errMsg
  use abortutils     , only : endrun
  use clm_varctl     , only : iulog
  use landunit_varcon, only : istice
  use clm_instur     , only : wt_lunit, wt_glc_mec
  use decompMod      , only : bounds_type, subgrid_level_gridcell, subgrid_level_column
  use filterColMod   , only : filter_col_type
  use ColumnType     , only : col
  use PatchType      , only : patch

  ! !PUBLIC TYPES:
  implicit none
  private
  save

  type, public :: glc_behavior_type
     private

     ! ------------------------------------------------------------------------
     ! Public data
     ! ------------------------------------------------------------------------

     ! If has_virtual_columns_grc(g) is true, then grid cell g has virtual columns for
     ! all possible glacier columns.
     !
     ! For the sake of coupling with CISM, this should only be needed within the icemask,
     ! where we need virtual columns for the sake of coupling with CISM. This is needed in
     ! order to (1) provide SMB in all elevation classes, in case it is being used with
     ! 1-way coupling (or to force a later TG run); (2) even with two-way coupling,
     ! provide SMB in the elevation classes above and below existing elevation classes,
     ! for the sake of vertical interpolation; (3) provide place-holder columns (which are
     ! already spun-up) for dynamic landunits; (4) ensure that all glacier columns are
     ! given spun-up initial conditions by init_interp.
     !
     ! More details on (4) (echoing the similar comment in subgridWeightsMod): We need all
     ! glacier and vegetated points to be active in the icemask region for the sake of
     ! init_interp - since we only interpolate onto active points, and we don't know which
     ! points will have non-zero area until after initialization (as long as we can't send
     ! information from glc to clm in initialization). (If we had an inactive glacier
     ! point in the icemask region, according to the weights on the surface dataset, and
     ! ran init_interp, this point would keep its cold start initialization values. Then,
     ! in the first time step of the run loop, it's possible that this point would become
     ! active because, according to glc, there is actually > 0% glacier in that grid
     ! cell. We don't do any state / flux adjustments in the first time step after
     ! init_interp due to glacier area changes, so this glacier column would remain at its
     ! cold start initialization values, which would be a Bad Thing. Ensuring that all
     ! glacier points within the icemask are active gets around this problem - as well as
     ! having other benefits, as noted above.)
     !
     ! However, by making this part of the user-modifiable "glc behavior", we make it easy
     ! for the user to add virtual columns, if this is desired for diagnostic
     ! purposes. One important reason why this may be desired is to produce coupler
     ! history forcings to force a later TG run, with SMB forcings outside the original
     ! CISM area. (Also, we cannot use icemask for all purposes, because it isn't known at
     ! initialization.)
     logical, allocatable, public :: has_virtual_columns_grc(:)

     ! If allow_multiple_columns_grc(g) is true, then grid cell g may have multiple
     ! glacier columns, for the different elevation classes. If
     ! allow_multiple_columns_grc(g) is false, then grid cell g is guaranteed to have at
     ! most one glacier column.
     logical, allocatable, public :: allow_multiple_columns_grc(:)

     ! If melt_replaced_by_ice_grc(g) is true, then any glacier ice melt in gridcell g
     ! runs off and is replaced by ice. Note that SMB cannot be computed in gridcell g if
     ! melt_replaced_by_ice_grc(g) is false, since we can't compute a sensible negative
     ! smb in that case.
     logical, allocatable, public :: melt_replaced_by_ice_grc(:)

     ! If ice_runoff_melted_grc(g) is true, then ice runoff generated by the
     ! CLM physics over glacier columns in gridcell g is melted (generating a negative
     ! sensible heat flux) and runs off as liquid. If it is false, then ice runoff is
     ! sent to the river model as ice (a crude parameterization of iceberg calving).
     logical, allocatable, public :: ice_runoff_melted_grc(:)

     ! ------------------------------------------------------------------------
     ! Private data
     ! ------------------------------------------------------------------------

     ! If collapse_to_atm_topo_grc(g) is true, then grid cell g has at most one glacier
     ! column, whose topographic height exactly matches the atmosphere's topographic
     ! height for that grid cell (so that there is no adjustment of atmospheric
     ! forcings).
     !
     ! Note that has_virtual_columns_grc(g) is guaranteed to be false if
     ! collapse_to_atm_topo_grc(g) is true.
     logical, allocatable :: collapse_to_atm_topo_grc(:)

   contains

     ! ------------------------------------------------------------------------
     ! Public routines
     ! ------------------------------------------------------------------------

     procedure, public  :: Init            ! version of Init meant for production use
     procedure, public  :: InitFromInputs  ! version of Init meant for unit testing (and called by other code in this class)
     procedure, public  :: InitSetDirectly ! version of Init meant for unit testing

     ! get number of subgrid units in glc landunit on one grid cell
     procedure, public  :: get_num_glc_subgrid

     ! returns true if memory should be allocated for the given glc column, and its
     ! weight on the landunit
     procedure, public  :: glc_col_exists

     ! returns true if glc columns on the given grid cell have dynamic type (type
     ! potentially changing at runtime)
     procedure, public  :: cols_have_dynamic_type

     ! Sets a column-level logical array to true for any ice column that needs
     ! downscaling, false for any ice column that does not need downscaling
     procedure, public  :: ice_cols_need_downscaling

     ! Sets a column-level logical array to true for any ice column that has
     ! dynamic type, false for any ice column that does not have dynamic type
     procedure, public :: cols_have_dynamic_type_array

     ! Sets a patch-level logical array to true for any ice column that has
     ! dynamic type, false for any ice column that does not have dynamic type
     procedure, public :: patches_have_dynamic_type_array

     ! update the column class types of any glc columns that need to be updated
     procedure, public  :: update_glc_classes

     ! ------------------------------------------------------------------------
     ! Public routines, for unit tests only
     ! ------------------------------------------------------------------------

     ! get the value of collapse_to_atm_topo at a given grid cell
     procedure, public :: get_collapse_to_atm_topo

     ! ------------------------------------------------------------------------
     ! Private routines
     ! ------------------------------------------------------------------------

     procedure, private :: InitAllocate

     ! reads GLACIER_REGION field from surface dataset
     procedure, private, nopass :: read_surface_dataset

     ! reads local namelist items
     procedure, private, nopass :: read_namelist

     ! returns a column-level filter of ice columns with the collapse_to_atm_topo
     ! behavior
     procedure, private :: collapse_to_atm_topo_ice_filterc

     ! update class of glc columns in regions where these are collapsed to a single
     ! column, given a filter
     procedure, private :: update_collapsed_columns_classes

  end type glc_behavior_type

  ! !PRIVATE MEMBER DATA:

  ! Longest name allowed for glacier_region_behavior, glacier_region_melt_behavior and
  ! glacier_region_ice_runoff_behavior
  integer, parameter :: max_behavior_name_len = 32

  ! Smallest and largest allowed values for a glacier region ID
  integer, parameter :: min_glacier_region_id = 0
  integer, parameter :: max_glacier_region_id = 10

  ! Value indicating that a namelist item has not been set
  character(len=*), parameter :: behavior_str_unset = 'UNSET'

  character(len=*), parameter, private :: sourcefile = &
       __FILE__

contains

  !-----------------------------------------------------------------------
  subroutine Init(this, begg, endg, NLFilename)
    !
    ! !DESCRIPTION:
    ! Initialize a glc_behavior_type object.
    !
    ! This version of Init is the one intended for production code use. It reads the
    ! information it needs from the surface dataset and namelist.
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    class(glc_behavior_type), intent(inout) :: this
    integer, intent(in) :: begg  ! beginning gridcell index
    integer, intent(in) :: endg  ! ending gridcell index
    character(len=*), intent(in)  :: NLFilename ! Namelist filename
    !
    ! !LOCAL VARIABLES:
    integer, allocatable :: glacier_region_map(:)
    character(len=max_behavior_name_len) :: glacier_region_behavior(min_glacier_region_id:max_glacier_region_id)
    character(len=max_behavior_name_len) :: glacier_region_melt_behavior(min_glacier_region_id:max_glacier_region_id)
    character(len=max_behavior_name_len) :: glacier_region_ice_runoff_behavior(min_glacier_region_id:max_glacier_region_id)

    character(len=*), parameter :: subname = 'Init'
    !-----------------------------------------------------------------------

    allocate(glacier_region_map(begg:endg))
    call this%read_surface_dataset(begg, endg, glacier_region_map(begg:endg))
    call this%read_namelist(NLFilename, glacier_region_behavior, &
         glacier_region_melt_behavior, glacier_region_ice_runoff_behavior)

    call this%InitFromInputs(begg, endg, &
         glacier_region_map(begg:endg), glacier_region_behavior, &
         glacier_region_melt_behavior, glacier_region_ice_runoff_behavior)

  end subroutine Init

  !-----------------------------------------------------------------------
  subroutine InitFromInputs(this, begg, endg, &
       glacier_region_map, glacier_region_behavior_str, glacier_region_melt_behavior_str, &
       glacier_region_ice_runoff_behavior_str)
    !
    ! !DESCRIPTION:
    ! Initialize a glc_behavior_type object given a map of glacier region IDs and an
    ! array of behavior specifications for each of these IDs.
    !
    ! This version should generally only be called directly by tests, but it is also used
    ! by the main production Init method.
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    class(glc_behavior_type), intent(inout) :: this
    integer, intent(in) :: begg  ! beginning gridcell index
    integer, intent(in) :: endg  ! ending gridcell index

    ! map of glacier region IDs
    integer, intent(in) :: glacier_region_map(begg:)

    ! string giving behavior for each glacier region ID
    ! allowed values are:
    ! - 'multiple': grid cells can potentially have multiple glacier elevation classes,
    !   but no virtual columns
    ! - 'virtual': grid cells have virtual columns: values are computed for every glacier
    !   elevation class, even those with 0 area
    ! - 'single_at_atm_topo': glacier landunits in these grid cells have a single column,
    !   whose elevation matches the atmosphere's topographic height (so that there is no
    !   adjustment due to downscaling)
    character(len=*), intent(in) :: glacier_region_behavior_str(min_glacier_region_id:)

    ! string giving treatment of ice melt for each glacier region ID
    ! allowed values are:
    ! - 'replaced_by_ice'
    ! - 'remains_in_place'
    character(len=*), intent(in) :: glacier_region_melt_behavior_str(min_glacier_region_id:)

    ! string giving treatment of ice runoff for each glacier region ID
    ! allowed values are:
    ! - 'remains_ice'
    ! - 'melted'
    character(len=*), intent(in) :: glacier_region_ice_runoff_behavior_str(min_glacier_region_id:)

    !
    ! !LOCAL VARIABLES:
    ! whether each glacier region ID is present in the glacier_region_map
    logical :: glacier_region_present(min_glacier_region_id:max_glacier_region_id)

    ! integer codes corresponding to glacier_region_behavior_str
    integer :: glacier_region_behavior(min_glacier_region_id:max_glacier_region_id)

    ! integer codes corresponding to glacier_region_melt_behavior_str
    integer :: glacier_region_melt_behavior(min_glacier_region_id:max_glacier_region_id)

    ! integer codes corresponding to glacier_region_ice_runoff_behavior_str
    integer :: glacier_region_ice_runoff_behavior(min_glacier_region_id:max_glacier_region_id)

    integer :: g
    integer :: my_id
    integer :: my_behavior
    integer :: my_melt_behavior
    integer :: my_ice_runoff_behavior

    ! possible glacier_region_behavior codes
    integer, parameter :: BEHAVIOR_MULTIPLE = 1
    integer, parameter :: BEHAVIOR_VIRTUAL = 2
    integer, parameter :: BEHAVIOR_SINGLE_AT_ATM_TOPO = 3

    ! possible glacier_region_melt_behavior codes
    integer, parameter :: MELT_BEHAVIOR_REPLACED_BY_ICE = 1
    integer, parameter :: MELT_BEHAVIOR_REMAINS_IN_PLACE = 2

    ! possible glacier_region_ice_runoff_behavior codes
    integer, parameter :: ICE_RUNOFF_BEHAVIOR_REMAINS_ICE = 1
    integer, parameter :: ICE_RUNOFF_BEHAVIOR_MELTED = 2

    ! value indicating that a behavior code has not been set (for glacier_region_behavior,
    ! glacier_region_melt_behavior or glacier_region_ice_runoff_behavior)
    integer, parameter :: BEHAVIOR_UNSET = -1

    character(len=*), parameter :: subname = 'InitFromInputs'
    !-----------------------------------------------------------------------

    SHR_ASSERT_ALL_FL((ubound(glacier_region_map) == (/endg/)), sourcefile, __LINE__)

    call check_glacier_region_map

    call determine_region_presence

    call translate_glacier_region_behavior
    call translate_glacier_region_melt_behavior
    call translate_glacier_region_ice_runoff_behavior

    call this%InitAllocate(begg, endg)

    do g = begg, endg
       my_id = glacier_region_map(g)
       my_behavior = glacier_region_behavior(my_id)
       my_melt_behavior = glacier_region_melt_behavior(my_id)
       my_ice_runoff_behavior = glacier_region_ice_runoff_behavior(my_id)

       ! This should only happen due to a programming error, not due to a user input error
       SHR_ASSERT_FL(my_behavior /= BEHAVIOR_UNSET, sourcefile, __LINE__)
       SHR_ASSERT_FL(my_melt_behavior /= BEHAVIOR_UNSET, sourcefile, __LINE__)
       SHR_ASSERT_FL(my_ice_runoff_behavior /= BEHAVIOR_UNSET, sourcefile, __LINE__)

       if (my_behavior == BEHAVIOR_VIRTUAL) then
          this%has_virtual_columns_grc(g) = .true.
       else
          this%has_virtual_columns_grc(g) = .false.
       end if

       if (my_melt_behavior == MELT_BEHAVIOR_REMAINS_IN_PLACE) then
          this%melt_replaced_by_ice_grc(g) = .false.
       else
          this%melt_replaced_by_ice_grc(g) = .true.
       end if

       if (my_ice_runoff_behavior == ICE_RUNOFF_BEHAVIOR_MELTED) then
          this%ice_runoff_melted_grc(g) = .true.
       else
          this%ice_runoff_melted_grc(g) = .false.
       end if

       ! For now, allow_multiple_columns_grc is simply the opposite of
       ! collapse_to_atm_topo_grc. However, we maintain the separate
       ! allow_multiple_columns_grc so that the public interface can stay the same if we
       ! differentiate between the two in the future - e.g., allowing for the possibility
       ! of a behavior where we have at most one glacier column, but not forced to the
       ! atmosphere's elevation.
       if (my_behavior == BEHAVIOR_SINGLE_AT_ATM_TOPO) then
          this%collapse_to_atm_topo_grc(g) = .true.
          this%allow_multiple_columns_grc(g) = .false.
       else
          this%collapse_to_atm_topo_grc(g) = .false.
          this%allow_multiple_columns_grc(g) = .true.
       end if
    end do

  contains
    subroutine check_glacier_region_map
      if (minval(glacier_region_map) < min_glacier_region_id) then
         write(iulog,*) subname//' ERROR: Expect GLACIER_REGION to be >= ', min_glacier_region_id
         write(iulog,*) 'minval = ', minval(glacier_region_map)
         call endrun(msg=' ERROR: GLACIER_REGION smaller than expected'// &
              errMsg(sourcefile, __LINE__))
      end if

      if (maxval(glacier_region_map) > max_glacier_region_id) then
         write(iulog,*) subname//' ERROR: Max GLACIER_REGION is ', &
              maxval(glacier_region_map)
         write(iulog,*) 'but max_glacier_region_id is only ', max_glacier_region_id
         write(iulog,*) 'Try increasing max_glacier_region_id in ', sourcefile
         call endrun(msg=' ERROR: GLACIER_REGION larger than expected'// &
              errMsg(sourcefile, __LINE__))
      end if
    end subroutine check_glacier_region_map

    subroutine determine_region_presence
      integer :: g
      integer :: my_id

      glacier_region_present(:) = .false.
      do g = begg, endg
         my_id = glacier_region_map(g)
         glacier_region_present(my_id) = .true.
      end do
    end subroutine determine_region_presence

    subroutine translate_glacier_region_behavior
      integer :: i

      do i = min_glacier_region_id, max_glacier_region_id
         glacier_region_behavior(i) = BEHAVIOR_UNSET

         if (glacier_region_present(i)) then
            SHR_ASSERT_ALL_FL((ubound(glacier_region_behavior_str) >= (/i/)), sourcefile, __LINE__)

            select case (glacier_region_behavior_str(i))
            case ('multiple')
               glacier_region_behavior(i) = BEHAVIOR_MULTIPLE
            case ('virtual')
               glacier_region_behavior(i) = BEHAVIOR_VIRTUAL
            case ('single_at_atm_topo')
               glacier_region_behavior(i) = BEHAVIOR_SINGLE_AT_ATM_TOPO
            case (behavior_str_unset)
               write(iulog,*) ' ERROR: glacier_region_behavior not specified for ID ', i
               write(iulog,*) 'You probably need to extend the glacier_region_behavior namelist array'
               call endrun(msg=' ERROR: glacier_region_behavior not specified for ID '// &
                    errMsg(sourcefile, __LINE__))
            case default
               write(iulog,*) ' ERROR: Unknown glacier_region_behavior for ID ', i
               write(iulog,*) glacier_region_behavior_str(i)
               write(iulog,*) 'Allowable values are: multiple, virtual, single_at_atm_topo'
               call endrun(msg=' ERROR: Unknown glacier_region_behavior'// &
                    errMsg(sourcefile, __LINE__))
            end select

         end if
      end do
    end subroutine translate_glacier_region_behavior

    subroutine translate_glacier_region_melt_behavior
      integer :: i

      do i = min_glacier_region_id, max_glacier_region_id
         glacier_region_melt_behavior(i) = BEHAVIOR_UNSET

         if (glacier_region_present(i)) then
            SHR_ASSERT_ALL_FL((ubound(glacier_region_melt_behavior_str) >= (/i/)), sourcefile, __LINE__)

            select case (glacier_region_melt_behavior_str(i))
            case ('replaced_by_ice')
               glacier_region_melt_behavior(i) = MELT_BEHAVIOR_REPLACED_BY_ICE
            case ('remains_in_place')
               glacier_region_melt_behavior(i) = MELT_BEHAVIOR_REMAINS_IN_PLACE
            case (behavior_str_unset)
               write(iulog,*) ' ERROR: glacier_region_melt_behavior not specified for ID ', i
               write(iulog,*) 'You probably need to extend the glacier_region_melt_behavior namelist array'
               call endrun(msg=' ERROR: glacier_region_melt_behavior not specified for ID '// &
                    errMsg(sourcefile, __LINE__))
            case default
               write(iulog,*) ' ERROR: Unknown glacier_region_melt_behavior for ID ', i
               write(iulog,*) glacier_region_melt_behavior_str(i)
               write(iulog,*) 'Allowable values are: replaced_by_ice, remains_in_place'
               call endrun(msg=' ERROR: Unknown glacier_region_melt_behavior'// &
                    errMsg(sourcefile, __LINE__))
            end select

         end if
      end do
    end subroutine translate_glacier_region_melt_behavior

    subroutine translate_glacier_region_ice_runoff_behavior
      integer :: i

      do i = min_glacier_region_id, max_glacier_region_id
         glacier_region_ice_runoff_behavior(i) = BEHAVIOR_UNSET

         if (glacier_region_present(i)) then
            SHR_ASSERT_ALL_FL((ubound(glacier_region_ice_runoff_behavior_str) >= (/i/)), sourcefile, __LINE__)

            select case (glacier_region_ice_runoff_behavior_str(i))
            case ('remains_ice')
               glacier_region_ice_runoff_behavior(i) = ICE_RUNOFF_BEHAVIOR_REMAINS_ICE
            case('melted')
               glacier_region_ice_runoff_behavior(i) = ICE_RUNOFF_BEHAVIOR_MELTED
            case (behavior_str_unset)
               write(iulog,*) ' ERROR: glacier_region_ice_runoff_behavior not specified for ID ', i
               write(iulog,*) 'You probably need to extend the glacier_region_ice_runoff_behavior namelist array'
               call endrun(msg=' ERROR: glacier_region_ice_runoff_behavior not specified for ID '// &
                    errMsg(sourcefile, __LINE__))
            case default
               write(iulog,*) ' ERROR: Unknown glacier_region_ice_runoff_behavior for ID ', i
               write(iulog,*) glacier_region_ice_runoff_behavior_str(i)
               write(iulog,*) 'Allowable values are: remains_ice, melted'
               call endrun(msg=' ERROR: Unknown glacier_region_ice_runoff_behavior'// &
                    errMsg(sourcefile, __LINE__))
            end select
         end if
      end do
    end subroutine translate_glacier_region_ice_runoff_behavior

  end subroutine InitFromInputs


  !-----------------------------------------------------------------------
  subroutine InitSetDirectly(this, begg, endg, &
       has_virtual_columns, collapse_to_atm_topo)
    !
    ! !DESCRIPTION:
    ! Initialize a glc_behavior_type object by directly setting has_virtual_columns and
    ! collapse_to_atm_topo
    !
    ! This version is meant for testing
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    class(glc_behavior_type), intent(inout) :: this
    integer, intent(in) :: begg  ! beginning gridcell index
    integer, intent(in) :: endg  ! ending gridcell index
    logical, intent(in) :: has_virtual_columns(begg:)
    logical, intent(in) :: collapse_to_atm_topo(begg:)
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'InitForTesting'
    !-----------------------------------------------------------------------

    SHR_ASSERT_ALL_FL((ubound(has_virtual_columns) == (/endg/)), sourcefile, __LINE__)
    SHR_ASSERT_ALL_FL((ubound(collapse_to_atm_topo) == (/endg/)), sourcefile, __LINE__)

    call this%InitAllocate(begg, endg)
    this%has_virtual_columns_grc(:) = has_virtual_columns(:)
    this%collapse_to_atm_topo_grc(:) = collapse_to_atm_topo(:)

  end subroutine InitSetDirectly


  !-----------------------------------------------------------------------
  subroutine InitAllocate(this, begg, endg)
    !
    ! !DESCRIPTION:
    ! Allocate variables in this object
    !
    ! !ARGUMENTS:
    class(glc_behavior_type), intent(inout) :: this
    integer, intent(in) :: begg  ! beginning gridcell index
    integer, intent(in) :: endg  ! ending gridcell index
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'InitAllocate'
    !-----------------------------------------------------------------------

    allocate(this%has_virtual_columns_grc (begg:endg)); this%has_virtual_columns_grc (:) = .false.
    allocate(this%allow_multiple_columns_grc(begg:endg)); this%allow_multiple_columns_grc(:) = .false.
    allocate(this%melt_replaced_by_ice_grc(begg:endg)); this%melt_replaced_by_ice_grc(:) = .false.
    allocate(this%collapse_to_atm_topo_grc(begg:endg)); this%collapse_to_atm_topo_grc(:) = .false.
    allocate(this%ice_runoff_melted_grc(begg:endg)); this%ice_runoff_melted_grc(:) = .false.

  end subroutine InitAllocate

  !-----------------------------------------------------------------------
  subroutine read_surface_dataset(begg, endg, glacier_region_map)
    !
    ! !DESCRIPTION:
    ! Reads GLACIER_REGION field from surface dataset, returns it in glacier_region_map
    !
    ! !USES:
    use clm_varctl , only : fsurdat
    use fileutils  , only : getfil
    use ncdio_pio  , only : file_desc_t, ncd_io, ncd_pio_openfile, ncd_pio_closefile
    use spmdMod    , only : masterproc
    use clm_varcon , only : grlnd
    !
    ! !ARGUMENTS:
    integer, intent(in)  :: begg  ! beginning grid cell index
    integer, intent(in)  :: endg  ! ending grid cell index
    integer, intent(out) :: glacier_region_map(begg:)
    !
    ! !LOCAL VARIABLES:
    integer, pointer :: glacier_region_map_ptr(:)  ! pointer version needed for ncd_io interface
    character(len=256) :: locfn        ! local filename
    type(file_desc_t)  :: ncid         ! netcdf id
    logical            :: readvar 

    character(len=*), parameter :: subname = 'read_surface_dataset'
    !-----------------------------------------------------------------------

    SHR_ASSERT_ALL_FL((ubound(glacier_region_map) == (/endg/)), sourcefile, __LINE__)

    if (masterproc) then
       write(iulog,*) 'Attempting to read GLACIER_REGION...'
    end if
    call getfil(fsurdat, locfn, 0)
    call ncd_pio_openfile(ncid, locfn, 0)
    allocate(glacier_region_map_ptr(begg:endg))
    call ncd_io(ncid=ncid, varname='GLACIER_REGION', flag='read', &
         data=glacier_region_map_ptr, dim1name=grlnd, readvar=readvar)
    if (.not. readvar) then
       call endrun(msg=' ERROR: GLACIER_REGION NOT on surfdata file'//errMsg(sourcefile, __LINE__))
    end if
    call ncd_pio_closefile(ncid)
    glacier_region_map(begg:endg) = glacier_region_map_ptr(begg:endg)
    deallocate(glacier_region_map_ptr)

  end subroutine read_surface_dataset

  !-----------------------------------------------------------------------
  subroutine read_namelist(NLFilename, glacier_region_behavior, &
       glacier_region_melt_behavior, glacier_region_ice_runoff_behavior)
    !
    ! !DESCRIPTION:
    ! Read local namelist items
    !
    ! !USES:
    use fileutils      , only : getavu, relavu, opnfil
    use shr_nl_mod     , only : shr_nl_find_group_name
    use clm_nlUtilsMod , only : find_nlgroup_name
    use spmdMod        , only : masterproc, mpicom
    use shr_mpi_mod    , only : shr_mpi_bcast
    !
    ! !ARGUMENTS:
    character(len=*), intent(in)  :: NLFilename ! Namelist filename
    character(len=max_behavior_name_len), intent(out) :: &
         glacier_region_behavior(min_glacier_region_id:max_glacier_region_id)
    character(len=max_behavior_name_len), intent(out) :: &
         glacier_region_melt_behavior(min_glacier_region_id:max_glacier_region_id)
    character(len=max_behavior_name_len), intent(out) :: &
         glacier_region_ice_runoff_behavior(min_glacier_region_id:max_glacier_region_id)
    !
    ! !LOCAL VARIABLES:
    integer :: unitn     ! unit for namelist file
    integer :: nml_error ! namelist i/o error flag

    character(len=*), parameter :: subname = 'read_namelist'
    !-----------------------------------------------------------------------

    namelist /clm_glacier_behavior/ &
         glacier_region_behavior, glacier_region_melt_behavior, &
         glacier_region_ice_runoff_behavior

    ! Initialize options to default values
    glacier_region_behavior(:) = behavior_str_unset
    glacier_region_melt_behavior(:) = behavior_str_unset
    glacier_region_ice_runoff_behavior(:) = behavior_str_unset

    if (masterproc) then
       unitn = getavu()
       call opnfil(NLFilename, unitn, 'F')
       call shr_nl_find_group_name(unitn, 'clm_glacier_behavior', status=nml_error)
       if (nml_error == 0) then
          read(unitn, nml=clm_glacier_behavior, iostat=nml_error)
          if (nml_error /= 0) then
             call endrun(msg='ERROR reading clm_glacier_behavior namelist'// &
                  errMsg(sourcefile, __LINE__))
          end if
       else
          call endrun(msg='ERROR finding clm_glacier_behavior namelist'// &
                      errMsg(sourcefile, __LINE__))
       end if
       call relavu( unitn )
    endif

    call shr_mpi_bcast(glacier_region_behavior, mpicom)
    call shr_mpi_bcast(glacier_region_melt_behavior, mpicom)
    call shr_mpi_bcast(glacier_region_ice_runoff_behavior, mpicom)

    if (masterproc) then
       write(iulog,*) ' '
       write(iulog,*) 'clm_glacier_behavior settings:'
       write(iulog,nml=clm_glacier_behavior)
       write(iulog,*) ' '
    end if

  end subroutine read_namelist


  !-----------------------------------------------------------------------
  subroutine get_num_glc_subgrid(this, gi, atm_topo, npatches, ncols, nlunits)
    !
    ! !DESCRIPTION:
    ! Get number of subgrid units in glc landunit on one grid cell
    !
    ! !USES:
    use clm_varpar      , only : maxpatch_glc
    !
    ! !ARGUMENTS:
    class(glc_behavior_type), intent(in) :: this
    integer , intent(in)  :: gi       ! grid cell index
    real(r8), intent(in)  :: atm_topo ! atmosphere's topographic height for this grid cell (m)
    integer , intent(out) :: npatches ! number of glacier_mec patches in this grid cell
    integer , intent(out) :: ncols    ! number of glacier_mec columns in this grid cell
    integer , intent(out) :: nlunits  ! number of glacier_mec landunits in this grid cell
    !
    ! !LOCAL VARIABLES:
    integer  :: m  ! loop index
    logical  :: col_exists
    real(r8) :: col_wt_lunit

    character(len=*), parameter :: subname = 'get_num_glc_subgrid'
    !-----------------------------------------------------------------------

    ncols = 0

    do m = 1, maxpatch_glc
       call this%glc_col_exists(gi = gi, elev_class = m, atm_topo = atm_topo, &
            exists = col_exists, col_wt_lunit = col_wt_lunit)
       if (col_exists) then
          ncols = ncols + 1
       end if
    end do

    if (this%collapse_to_atm_topo_grc(gi) .and. &
         wt_lunit(gi, istice) > 0.0_r8) then
       ! For grid cells with the collapse_to_atm_topo behavior, with a non-zero weight
       ! ice landunit, we expect exactly one column
       SHR_ASSERT_FL(ncols == 1, sourcefile, __LINE__)
    end if

    if (ncols > 0) then
       npatches = ncols
       nlunits = 1
    else
       npatches = 0
       nlunits = 0
    end if
  
  end subroutine get_num_glc_subgrid

  !-----------------------------------------------------------------------
  subroutine glc_col_exists(this, gi, elev_class, atm_topo, exists, col_wt_lunit)
    !
    ! !DESCRIPTION:
    ! For the given glc column, with elevation class index elev_class, in grid cell
    ! gi: sets exists to true if memory should be allocated for this column, and sets
    ! col_wt_lunit to the column's weight on the ice landunit.
    !
    ! If exists is false, then col_wt_lunit is arbitrary and should be ignored.
    !
    ! !USES:
    use glc_elevclass_mod, only : glc_get_elevation_class, GLC_ELEVCLASS_ERR_NONE
    use glc_elevclass_mod, only : GLC_ELEVCLASS_ERR_TOO_LOW, GLC_ELEVCLASS_ERR_TOO_HIGH
    use glc_elevclass_mod, only : glc_errcode_to_string
    !
    ! !ARGUMENTS:
    class(glc_behavior_type), intent(in) :: this
    integer,  intent(in)  :: gi           ! grid cell index
    integer,  intent(in)  :: elev_class   ! elevation class index
    real(r8), intent(in)  :: atm_topo     ! atmosphere's topographic height for this grid cell (m)
    logical,  intent(out) :: exists       ! whether memory should be allocated for this column
    real(r8), intent(out) :: col_wt_lunit ! column's weight on the ice landunit
    !
    ! !LOCAL VARIABLES:
    integer :: atm_elev_class ! elevation class corresponding to atmosphere topographic height
    integer :: err_code

    character(len=*), parameter :: subname = 'glc_col_exists'
    !-----------------------------------------------------------------------

    ! Set default outputs
    exists = .false.
    col_wt_lunit = wt_glc_mec(gi, elev_class)

    if (this%collapse_to_atm_topo_grc(gi)) then
       if (wt_lunit(gi, istice) > 0.0_r8) then
          call glc_get_elevation_class(atm_topo, atm_elev_class, err_code)
          if ( err_code == GLC_ELEVCLASS_ERR_NONE .or. &
               err_code == GLC_ELEVCLASS_ERR_TOO_LOW .or. &
               err_code == GLC_ELEVCLASS_ERR_TOO_HIGH) then
             ! These are all acceptable "errors" - it is even okay for these purposes if
             ! the elevation is lower than the lower bound of elevation class 1, or
             ! higher than the upper bound of the top elevation class.

             ! Do nothing
          else
             write(iulog,*) subname, ': ERROR getting elevation class for topo = ', atm_topo
             write(iulog,*) glc_errcode_to_string(err_code)
             call endrun(subgrid_index=gi, subgrid_level=subgrid_level_gridcell, &
                  msg=subname//': ERROR getting elevation class')
          end if

          if (elev_class == atm_elev_class) then
             exists = .true.
             col_wt_lunit = 1._r8
          else
             exists = .false.
             col_wt_lunit = 0._r8
          end if
       end if

    else  ! collapse_to_atm_topo_grc .false.
       if (this%has_virtual_columns_grc(gi)) then
          exists = .true.
       else if (wt_lunit(gi, istice) > 0.0_r8 .and. &
            wt_glc_mec(gi, elev_class) > 0.0_r8) then
          ! If the landunit has non-zero weight on the grid cell, and this column has
          ! non-zero weight on the landunit...
          exists = .true.
       end if
    end if

  end subroutine glc_col_exists

  !-----------------------------------------------------------------------
  function cols_have_dynamic_type(this, gi)
    !
    ! !DESCRIPTION:
    ! Returns true if glc columns on the given grid cell have dynamic type (i.e.,
    ! type potentially changing at runtime)
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    logical :: cols_have_dynamic_type  ! function result
    class(glc_behavior_type), intent(in) :: this
    integer, intent(in) :: gi  ! grid cell index
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'cols_have_dynamic_type'
    !-----------------------------------------------------------------------

    if (this%collapse_to_atm_topo_grc(gi)) then
       cols_have_dynamic_type = .true.
    else
       cols_have_dynamic_type = .false.
    end if

  end function cols_have_dynamic_type

  !-----------------------------------------------------------------------
  subroutine ice_cols_need_downscaling(this, bounds, num_icec, filter_icec, &
       needs_downscaling_col)
    !
    ! !DESCRIPTION:
    ! Sets needs_downscaling_col to true for any ice column that needs downscaling,
    ! false for any ice column that does not need downscaling.
    !
    ! Outside of filter_icec, leaves needs_downscaling_col untouched.
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    class(glc_behavior_type) , intent(in) :: this
    type(bounds_type)        , intent(in) :: bounds
    integer                  , intent(in) :: num_icec       ! number of points in filter_icec
    integer                  , intent(in) :: filter_icec(:) ! col filter for ice
    logical                  , intent(inout) :: needs_downscaling_col( bounds%begc: )
    !
    ! !LOCAL VARIABLES:
    integer :: fc
    integer :: c
    integer :: g

    character(len=*), parameter :: subname = 'ice_cols_need_downscaling'
    !-----------------------------------------------------------------------

    SHR_ASSERT_ALL_FL((ubound(needs_downscaling_col) == (/bounds%endc/)), sourcefile, __LINE__)

    do fc = 1, num_icec
       c = filter_icec(fc)
       g = col%gridcell(c)

       if (this%collapse_to_atm_topo_grc(g)) then
          needs_downscaling_col(c) = .false.
       else
          needs_downscaling_col(c) = .true.
       end if
    end do

  end subroutine ice_cols_need_downscaling

  !-----------------------------------------------------------------------
  subroutine cols_have_dynamic_type_array(this, begc, endc, has_dynamic_type_col)
    !
    ! !DESCRIPTION:
    ! Sets a column-level logical array to true for any ice column that has
    ! dynamic type, false for any ice column that does not have dynamic type.
    !
    ! The value is undefined for non-ice columns.
    !
    ! !ARGUMENTS:
    class(glc_behavior_type) , intent(in) :: this
    integer, intent(in) :: begc  ! beginning column index
    integer, intent(in) :: endc  ! ending column index
    logical, intent(inout) :: has_dynamic_type_col( begc: )
    !
    ! !LOCAL VARIABLES:
    integer :: c
    integer :: g

    character(len=*), parameter :: subname = 'cols_have_dynamic_type_array'
    !-----------------------------------------------------------------------

    SHR_ASSERT_ALL_FL((ubound(has_dynamic_type_col) == (/endc/)), sourcefile, __LINE__)

    do c = begc, endc
       g = col%gridcell(c)
       ! Users shouldn't rely on the values set for non-ice columns, but it's simpler
       ! just to set this for all column types.
       if (this%collapse_to_atm_topo_grc(g)) then
          has_dynamic_type_col(c) = .true.
       else
          has_dynamic_type_col(c) = .false.
       end if
    end do

  end subroutine cols_have_dynamic_type_array

  !-----------------------------------------------------------------------
  subroutine patches_have_dynamic_type_array(this, begp, endp, has_dynamic_type_patch)
    !
    ! !DESCRIPTION:
    ! Sets a patch-level logical array to true for any ice patch that has
    ! dynamic type, false for any ice patch that does not have dynamic type.
    !
    ! The value is undefined for non-ice patches.
    !
    ! !ARGUMENTS:
    class(glc_behavior_type) , intent(in) :: this
    integer, intent(in) :: begp  ! beginning patch index
    integer, intent(in) :: endp  ! ending patch index
    logical, intent(inout) :: has_dynamic_type_patch( begp: )
    !
    ! !LOCAL VARIABLES:
    integer :: p
    integer :: g

    character(len=*), parameter :: subname = 'patches_have_dynamic_type_array'
    !-----------------------------------------------------------------------

    SHR_ASSERT_ALL_FL((ubound(has_dynamic_type_patch) == (/endp/)), sourcefile, __LINE__)

    do p = begp, endp
       g = patch%gridcell(p)
       ! Users shouldn't rely on the values set for non-ice patches, but it's simpler
       ! just to set this for all patch types.
       if (this%collapse_to_atm_topo_grc(g)) then
          has_dynamic_type_patch(p) = .true.
       else
          has_dynamic_type_patch(p) = .false.
       end if
    end do

  end subroutine patches_have_dynamic_type_array

  !-----------------------------------------------------------------------
  subroutine update_glc_classes(this, bounds, topo_col)
    !
    ! !DESCRIPTION:
    ! Update the column class types of any glc columns that need to be updated.
    !
    ! Assumes that topo_col has already been set appropriately.
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    class(glc_behavior_type), intent(in) :: this
    type(bounds_type), intent(in) :: bounds
    real(r8), intent(in) :: topo_col( bounds%begc: )
    !
    ! !LOCAL VARIABLES:
    type(filter_col_type) :: collapse_filterc

    character(len=*), parameter :: subname = 'update_glc_classes'
    !-----------------------------------------------------------------------

    collapse_filterc = this%collapse_to_atm_topo_ice_filterc(bounds)
    call this%update_collapsed_columns_classes(bounds, collapse_filterc, topo_col)

  end subroutine update_glc_classes

  !-----------------------------------------------------------------------
  subroutine update_collapsed_columns_classes(this, bounds, collapse_filterc, topo_col)
    !
    ! !DESCRIPTION:
    ! Update class of glc columns in regions where these are collapsed to a single
    ! column, given a filter.
    !
    ! Assumes that topo_col has already been updated appropriately for these columns.
    !
    ! !USES:
    use glc_elevclass_mod, only : glc_get_elevation_class, GLC_ELEVCLASS_ERR_NONE
    use glc_elevclass_mod, only : GLC_ELEVCLASS_ERR_TOO_LOW, GLC_ELEVCLASS_ERR_TOO_HIGH
    use glc_elevclass_mod, only : glc_errcode_to_string
    use column_varcon    , only : ice_class_to_col_itype
    !
    ! !ARGUMENTS:
    class(glc_behavior_type), intent(in) :: this
    type(bounds_type), intent(in) :: bounds
    type(filter_col_type), intent(in) :: collapse_filterc
    real(r8), intent(in) :: topo_col( bounds%begc: )
    !
    ! !LOCAL VARIABLES:
    integer :: fc         ! filter index
    integer :: c          ! column index
    integer :: elev_class ! elevation class of the single column on the ice landunit
    integer :: err_code
    integer :: new_coltype

    character(len=*), parameter :: subname = 'update_collapsed_columns_classes'
    !-----------------------------------------------------------------------

    SHR_ASSERT_ALL_FL((ubound(topo_col) == (/bounds%endc/)), sourcefile, __LINE__)

    do fc = 1, collapse_filterc%num
       c = collapse_filterc%indices(fc)

       call glc_get_elevation_class(topo_col(c), elev_class, err_code)
       if ( err_code == GLC_ELEVCLASS_ERR_NONE .or. &
            err_code == GLC_ELEVCLASS_ERR_TOO_LOW .or. &
            err_code == GLC_ELEVCLASS_ERR_TOO_HIGH) then
          ! These are all acceptable "errors" - it is even okay for these purposes if
          ! the elevation is lower than the lower bound of elevation class 1, or
          ! higher than the upper bound of the top elevation class.
          
          ! Do nothing
       else
          write(iulog,*) subname, ': ERROR getting elevation class for topo = ', &
               topo_col(c)
          write(iulog,*) glc_errcode_to_string(err_code)
          call endrun(subgrid_index=c, subgrid_level=subgrid_level_column, &
               msg=subname//': ERROR getting elevation class')
       end if

       new_coltype = ice_class_to_col_itype(elev_class)
       if (new_coltype /= col%itype(c)) then
          call col%update_itype(c = c, itype = new_coltype)
       end if
    end do
          
  end subroutine update_collapsed_columns_classes

  !-----------------------------------------------------------------------
  function collapse_to_atm_topo_ice_filterc(this, bounds) result(filter)
    !
    ! !DESCRIPTION:
    ! Returns a column-level filter of ice columns with the collapse_to_atm_topo behavior
    !
    ! !USES:
    use filterColMod, only : filter_col_type, col_filter_from_grcflags_ltypes
    !
    ! !ARGUMENTS:
    class(glc_behavior_type), intent(in) :: this
    type(filter_col_type) :: filter  ! function result
    type(bounds_type), intent(in) :: bounds
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'collapse_to_atm_topo_ice_filterc'
    !-----------------------------------------------------------------------

    ! Currently this creates the filter on the fly, recreating it every time this
    ! function is called. Since this is a static filter, we could just compute it once
    ! and save it, returning the already-computed filter when this function is called.
    ! However, the problem with that is the need to have a different filter for each
    ! clump (and potentially another filter for calls from outside a clump loop). This
    ! will become easier to handle if we rework CLM's threading so that there is a
    ! separate instance of each object for each clump: in that case, we'll have multiple
    ! instances of glc_behavior_type, each corresponding to one clump, each with its own
    ! filter.

    filter = col_filter_from_grcflags_ltypes( &
         bounds = bounds, &
         grcflags = this%collapse_to_atm_topo_grc(bounds%begg:bounds%endg), &
         ltypes = [istice], &
         include_inactive = .true.)

  end function collapse_to_atm_topo_ice_filterc

  !-----------------------------------------------------------------------
  function get_collapse_to_atm_topo(this, gi) result(collapse_to_atm_topo)
    !
    ! !DESCRIPTION:
    ! Get the value of collapse_to_atm_topo at a given grid cell
    !
    ! This function just exists to support unit testing, and should not be called from
    ! production code.
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    logical :: collapse_to_atm_topo  ! function result
    class(glc_behavior_type), intent(in) :: this
    integer, intent(in) :: gi ! grid cell index
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'get_collapse_to_atm_topo'
    !-----------------------------------------------------------------------

    collapse_to_atm_topo = this%collapse_to_atm_topo_grc(gi)

  end function get_collapse_to_atm_topo

end module glcBehaviorMod
